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Abstract 

A formula is given for the propagation of errors during matrix in- 
version. An explicit calculation for a 2 x 2 matrix using both the 
formula and a Monte Carlo calculation are compared. A prescrip- 
tion is given to determine when a matrix with uncertain elements is 
sufficiently nonsingular for the calculation of the covariances of the 
inverted matrix elements to be reliable. 
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1 Introduction 

There are many problems that involve solving a set of simultaneous linear 
equations. In many instances, the coefficients of the variables are uncertain 
and these uncertainties need to be taken into account in the final solution. 
The results presented in this paper are general and can be applied to any 
problem involving the solution of linear equations. Moreover the full covari- 
ance matrix is calculated. The use of the covariance formula is illustrated 
using an example where a number of branching ratios of a given particle were 
simultaneously measured. 

A recent paper reporting results on the decay of the tau lepton (|I]) used 
a matrix inversion technique to solve for several branching ratios simultane- 
ously. In that paper, the treatment of the statistical uncertainties assumed 
only diagonal errors. This paper develops a formula for the covariance of 
the inverse matrix elements. The treatment follows the propagation of errors 
formalism (Q) for small errors. A full treatment of the covariance matrix is 
necessary because the off-diagonal errors can be significant. 
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The fraction of time a particle decays into a given final state is defined to 
be the branching ratio. The branching ratios of a particle decaying into a set 
of m possible final states may be determined by finding m selection criteria 
that have an efficiency for selecting the desired decay channel. The selection 
criteria are usually not fully efficient so a set of linear equations results, 

eii-Bi + £12-82 + £13-83 + ■ ■ ■ + eim-Bm = fi 

621-81 + €22-82 + £23-83 + ■ ■ ■ + t2mBm = /2 

£31-81 + £32-83 + £33-83 + ■ ■ ■ + £3m-8m = fs (1) 
£ml-8l + £m2-82 + ^m^B^ + ■ ■ ■ + emmBm = fm 

where Bj is the j^^ unknown branching ratio, fi is the fraction of events 
chosen by the i^^ selection from a measured sample with an efficiency e^j. The 
fractions are usually corrected for background using Monte Carlo estimates. 

The efficiency is typically calculated by applying the same selections to 
a Monte Carlo sample of events. Thus, the efficiency is given by 

where n^*^ is the number of events of decay channel j selected by selection 
i and A'j^^ is the total number events of type j in the Monte Carlo sample. 
As can be seen from equation |^, the elements of e^j are positive definite and 
less than or equal to one. An uncertainty for each efficiency matrix element 
can be determined from the Monte Carlo statistics. 

The simultaneous equations given by equation |l| can be written in matrix 
form, 

eB = f. (3) 
The vector of branching ratios can be solved from the matrix equation, 

B = e-'f , (4) 

where is the inverse of the efficiency matrix. The matrix e must be 
nonsingular in order for the inverse to exist or equivalently the determinant 
must be nonzero. A singular e signals an ill defined set of selections or event 
types. 



2 



The uncertainty of the branching ratios, B, can be written in terms of the 
uncertainties on the elements of the matrices / and e~^. The errors for / are 
determined from the data and the estimates of the background. The errors 
on can be calculated in terms of the covariances of the elements of e. The 
error formulae for used in the literature are described in section 2. A 
more general treatment of the errors taking into account the full covariance 
matrix is presented in section 3 whilst the detailed derivations are given in 
Appendix A. 

An analytical example based on matrices of order 2 x 2 is developed in 
Appendix B and a Monte Carlo study based on the 2x2 case is discussed 
in sections 4 and 5. The Monte Carlo studies will illustrate that the formula 
is only realistic when the efficiency matrix is sufficiently far away from being 
singular. The calculation of the covariances becomes more reliable as the 
ratio of the value of the determinant of the efficiency matrix, to the value 
of the uncertainty on the determinant becomes larger. A formula for the 
uncertainty of a determinant is derived in Appendix C. The paper concludes 
with some general observations. 

2 Calculation of Errors 

Errors are often estimated by ignoring the off-diagonal elements of the covari- 
ance matrix. This is the correct procedure if the quantities are independent 
of each other. A calculation of the uncertainty on the branching ratios, as- 
suming the elements of the inverted matrix are statistically independent of 
each other, yields Q 

where the uncertainties are denoted by a and a sum over repeated indices is 
assumed unless otherwise noted. The errors on the elements of the inverted 
efficiency matrix, a^-i, are calculated from the known errors on the efficiency 
matrix. 

The uncertainties on the inverse efficiencies, have been determined (see 
for example reference (|TD) by differentiating the matrix equation ee~^ = / 
and then applying a matrix analogy to the error propagation formula to yield, 

[cr,-i]^j- =1 [e~"^]jm[0"e]mn[e~"^]nj ^ ) (6) 

^Square brackets have been used to separate the superscripts and subscripts used to 
identify the quantity of interest from any indices and mathematical operations. 
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which has the same form as a transformation to a new basis (^. Equa- 
tion P is a reasonable approximation but is not correct. It also neglects any 
correlations between the elements of the inverse matrix, e~^, which can be 
significant, as shown in the next section. 



3 The Covariance of the Elements of e ^ 

Matrix inversion is a nonlinear operation. It is always possible to write 
the inverse of a matrix in terms of the matrix of cofactors divided by the 
determinant (|^). One sees explicitly in Appendix C that each element of an 
inverse matrix has elements of the original matrix in common. Therefore the 
inverse matrix elements clearly are correlated. 

Consider the matrix elements e^j with uncertainties given in the most 
general case by the covariances cov(eQ/3, eab)- The inverse matrix elements 
e^-^, in general, have covariances cov(e~^, e~j,^), which can be written as, 

cov(e-^, e-,^) = e-lejje-^eif^^cov{eij, eki) ■ (7) 

The full derivation of this equation is given in Appendix A (see equation ^ 
plOl) . The usual case where there are no correlations between the elements of 
the efficiency matrix (see equations |A-20| , |A-21| and [A-22| ) is given by 



cov {tij^tki) = [(r^]^ij6ikSji (no summation). (8) 
Hence the full set of covariances of are given by 

cov(e-^, e-J) = ([e-^.[e-']a.)K]J([e~^],,[e-^-.) , (9) 

where there is no sum in this case over repeated indices inside the parentheses. 
The variance of an element of the inverse efficiency matrix can be written as, 

[a,^.]l, ^ covie-l^e-J,) = [e-']l,K]l[e-% . (10) 

This equation is exact and replaces the approximate equation ^. Note that 



each term of equation |10| is squared before making the sum whereas in equa- 
tion ^ the sum is done first. 

The complete expression for the uncertainties on the branching ratios can 
then be calculated from, 

cov{Bi, Bj) = /«//3Cov(erJ, ej^) + er^hji^cov{fk, fi) , (11) 
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where the measured fractions of events are often uncorrelated, 

cov{fkJi) = [(Jf]lSM (no summation). (12) 
Equation ^ is a generahzation of equation ^, including all of the covariances. 

4 A Monte Carlo Study - Part 1 

The properties of equation || have been studied using a Monte Carlo simula- 
tion of a 2 X 2 matrix. Relatively small errors have been used to satisfy the 
small error approximation used to derive the error propagation formula (^. 
The results of the Monte Carlo calculations directly can be compared to the 
analytic formulae developed in Appendix B. 

Given an initial efficiency matrix e, and the variance on each element, 
[celfj, the m*'* random instance of the matrix is generated by 

where is a normally distributed pseudorandom deviate with a mean of 
zero and a standard deviation of one. A set of matrices are created, are 
then inverted. The covariances are calculated from 

cov(e;^, = {e-l e-,') - (e-^) {e-,') , (14) 

where ( ) is the mean of the quantity enclosed in the brackets. 
A sample of A^ = 10, 000 instances of the matrix e± where, 

_ / 0.700 0.200 \ f 0.007 0.002 \ , . 

^ ~ 0.400 0.600 ) \ 0.004 0.006 j ' ^ ^ 

were generated with each element of e normally distributed about the central 
values with the one percent errors, a^, indicated. The inverse of the matrix 
e is given by 

_i _ / 1.765 -0.588 
^ ~ 1^ -1.177 2.059 

and the determinant, dete = 0.340 The covariances can be calculated using 
only the inverse matrix, e~^, and the error matrix, a^. The numerical values 
for the 16 covariances of e~^, calculated from the Monte Carlo and analyti- 
cally using equation ^, are given in Table |I|. The symmetry of the covariance 
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-2.642 


2.514 


-2.619 


-0.007 


0.009 
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Table 1: The covariance matrix for each element of calculated by the 
Monte Carlo and calculated using the analytic formula. The last column is 
the fractional difference between the Monte Carlo and analytic calculation. 



operation implies there should be 10 unique numbers as is apparent from the 
table. 

Figure |1| shows a contour plot of the element ej~2^ versus ej~/ generated 
from the Monte Carlo. The contours are approximate lines of constant prob- 
abihty density. The entries are distributed around the values of the calcu- 
lated inverse matrix elements and are clearly correlated. The slope of the 
error ellipses has been calculated by substituting the analytical values of the 
covariances from Table |^ into a formula derived from reference (|^) . The cal- 
culated slope is in good agreement with the slope of the major axis of the 
error ellipses. 

Recall that the efficiency matrix must be nonsingular in order to invert 
it. In the Monte Carlo calculation the individual elements of e are varied and 
hence, if the variations are large enough, it is possible for the determinant 
to become zero. Figure ^ shows a histogram of the Monte Carlo calculation 
of the determinant. The mean value is 0.340 as expected with a root mean 
square deviation of 0.006 . The shape of the distribution is normal, see 
figure 1^, also as expected given the determinant is the sum of products of 
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Figure 1: The 10,000 elements ej"2 and ef^ generated in tiie Monte Carlo 
simulation form a two dimensional distribution shown here as a contour plot. 
The contours are at the level of 50, 100, 150, 200 and 250 counts per bin (each 
bin is 6.25 x 10^^ x 6.25 x 10^'^). The correlation between the two inverse 
matrix elements is clear. The data are centred around the exact calculation 
of the inverse, shown by the large "+" sign. The slope of the line is the 
theoretical calculation of the correlation coefficient. It is in good agreement 
with the slope of the axis of the approximately elliptical contours. 
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normally distributed quantities. In this particular case the mean value of 
the determinant is 56.7 standard deviations from zero. Increasing the errors 
or reducing the determinant will both produce a higher likelihood of the 
determinant fluctuating nearer to zero. 

As noted earlier the inverse matrix elements all have a common factor 
of l/|e|. Therefore values of the determinant near zero produce very large 
values for the inverse matrix elements. Moreover, the reciprocal of a normal 
distribution is not normally distributed but is instead asymmetric with a tail 
towards large values. In figure one sees that any tails are insignificant for 
determinants many standard deviations away from zero. 

Figure |^ shows an identical analysis for the matrix, 

, _ / 0.400 0.500 \ f 0.004 0.005 \ . . 

^ ~ 0.400 0.600 J \ 0.004 0.006 ) ' ^ ' 

In this case det e' = 0.040 and its root mean square deviation, due to the one 
percent errors, is 0.004 . Now the distribution of the reciprocal determinants, 
seen in Fugure is clearly skewed. Table ^ lists the values of the covariance 
matrix calculated for each element of the inverse matrix, e'~^, using the Monte 



Carlo method, equation and the analytic expression, equation ||. 

In Table ^ the percent difference between the Monte Carlo and the ana- 
lytic formula has a mean value of approximately one percent averaged over 
the 10 independent values of the covariance. This is consistent with the pre- 
cision expected for a 10,000 event Monte Carlo and demonstrates that the 
analytic expression reproduces the covariances well for the far from singular 
case. However, the same mean for the values in Table |^ is eleven percent. 
Moreover, the magnitude of the Monte Carlo value is larger than the analytic 
value for each covariance. The analytic expression becomes less accurate and 
no longer appropriate to use as the matrix becomes closer to being singular. 



5 A Monte Carlo Study - Part 2 

In order to explore the range of uncertainty on how well the analytic expres- 
sion models the covariances, the Monte Carlo simulation described in Part 1 
was modified. Instead of specifying a particular 2x2 matrix for analysis, 5000 
2x2 matrices were generated by choosing the elements of each matrix from 
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Figure 2: The histogram labelled (a) shows the determinant of the matrix e 
described in the text. The values of the determinant are clearly not near zero 
and consequently the histogram of the reciprocal of the determinant shown 
in (b) is well behaved. 
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Figure 3: The histogram labelled (a) shows the determinant of the matrix 
e described in the text. The values of the determinant approach zero and 
consequently an asymmetry is observed in the histogram, labelled (b) , of the 
reciprocal distribution. 
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Covariance Matrix e' ^ 


Monte Carlo x(10-^) 


Analytic x(10"^) 


Fractional Difference 




2.812 -2.554 
-2.040 1.855 


2.498 -2.269 
-1.815 1.650 


0.112 0.112 
0.110 0.110 




-2.554 2.337 
1.855 -1.699 


-2.269 2.078 
1.650 -1.513 


0.112 0.111 
0.110 0.110 


cov(e2/,e-/) 


-2.040 1.855 
1.492 -1.358 


-1.815 1.650 
1.330 -1.210 


0.110 0.110 
0.109 0.109 


cov(e22\e-/) 


1.855 -1.699 
-1.358 1.244 


1.650 -1.513 
-1.210 1.110 


0.110 0.111 
0.109 0.108 



Table 2: The covariance matrix for each element of e'^^ calculated by the 
Monte Carlo and calculated using the analytic formula. The last column is 
the fractional difference between the Monte Carlo and analytic calculation. 



a uniform distribution in the interval [0, 1] . Each of the randomly gener- 
ated matrices was then analyzed using exactly the same procedure previously 
discussed. One percent errors were used as before. 

None of the matrix inversions actually failed. However, many cases in- 
volved determinants close to zero. Figure ^(a) shows the fractional difference 
between the Monte Carlo and the analytic calculation for cov(en^, £21^), 

COv(ef/, e2i"^)MC - COv{eil, e2i^)Analytic ON 

cov(eii , €21 )mc 

plotted against the determinant of e divided by its uncertainty, dete/(T|e|. 
The number of sigmas the determinant may be from zero can be shown to be 
±70.7 for one percent errors. The fractional differences between the Monte 
Carlo and analytic evaluations are distributed around approximately zero for 
large numbers of sigma. However, near 10 sigmas the Monte Carlo calculated 
covariance becomes very large and the distribution becomes centred around 
one. Under the scatter plot is figure |5(b). This is a plot of the mean value of 
the fractional difference plotted against the number of sigmas the determi- 
nant is away from zero. The error bars are the root mean square deviation of 
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the distribution. The means for —10. < dete/(T|e| < 10. are not plotted. The 
largest deviation, for the bin 10-20, is less than four percent. The sigmas are 
approximately two percent. Hence, the analytic formula does a very good 
job of estimating the covariance as long as the determinant of the matrix is 
at least 10 sigmas from zero. A simple expression to calculate the error on a 
determinant is given in Appendix C. 

6 Conclusions 

The covariances of inverse matrix elements are nonzero in general. There- 
fore, the propagation of errors for formulae that depend on inverted matrices 
requires using the covariances of the inverse matrix elements. A concise for- 
mula is developed in the small error limit for the covariance of the inverted 
matrix elements. It is shown to work for nonsingular matrices. In particular 
the determinant must not become zero when the matrix elements are allowed 
to vary within their uncertainties. 

An attempt was made to study how well the formula estimates the co- 
variances for a random set of 2 x 2 matrices with positive elements between 
zero and one. On the average the formula was accurate to better than four 
percent with a root mean square deviation equal to two percent for matrices 
that have determinants more than ten sigmas from zero. 

7 Acknowledgements 

This project was supported by grants from the Natural Science and Engi- 
neering Research Council of Canada. 

References 

[1] The OPAL Collaboration, K. Ackerstaff et al., Eur. Phys. J. 
€4(1998)193-206. 

[2] The Particle Data Group, C. Caso et al.. Section 29, Eur. Phys. J. 
C3(1998) 1-794. 

[3] L. Lyons, Statistics for Nuclear and Particle Physics, Cambridge Univer- 
sity Press, 1989. 



12 



1 


1 1 1 1 


(a) 


u 

2 0.75 






:tional 

In 






2 0.25 











-60 -40 -20 20 40 60 

det(8)/o„ 



-0.1 



1 1 1 r 




^ (b) 



J I I L 



-60 -40 -20 20 40 60 

det(E)/a_ 
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calculations as a function of det e/a\e\ for cov(ej"/, €21) is plotted in part (a). 
The fractional difference is normalized by the Monte Carlo value. In part (b) 
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slices of det e/a\e\ between -70 and 70. The two centre shces are not plotted. 
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A Derivation of the Covariance 

The formal solution to the system of real linear equations, 

Yi^AijXj i,j^l...,n (A-1) 

where are the unknowns and A is a real n x n matrix, is given by 

X, = AT^Yj . (A-2) 

A sum over repeated indices is assumed unless otherwise noted. The uncer- 
tainties on the elements of A are assumed to be known and given by 

C0Y{Aij,Aaf3) , (A-3) 

the covariance of Aij and Aq,^. If each element of A is an independent random 
quantity, then only the terms cov{Aij, Aij) are nonzero. In other words, each 
element of the matrix will have an associated variance. 

Similarly the uncertainties on the elements of Y are assumed to be known 
and given by cov{Yi,Yj). Again, if the elements of Y are independent, the 
off-diagonal terms will be zero. 

The most general form of the covariance for the Xi is given by the error 
propagation formula, 

^ dXj dXi / . 1 . OXj QXa ,,,,,, , . 
cov(X., X,) = ^^cov(^;j, A-J) + oov(n, Y,) . (A-4) 

The partial derivatives are given by: 

dXi 



dXj 

dYk 
dXj 



dYi 



SaiYp , (A-5) 

SajY, , (A-6) 

A^' , (A-7) 

Aji' , (A-8) 
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where Sij is the Kronecker delta. Therefore, after relabelhng summed-over 
indices, equation \A-4\ becomes 

cov(X„ X,) = y;F^cov(4ri, AJ^) + A-j}Aj,\oY{Yk, Yi) . (A-9) 

The only unknown quantity is cov {A^^ , Aj^) . 

The inverse matrix elements can be considered as functions of the original 
matrix elements, 

A-'p = A-JiA,) . (A-10) 
The error propagation formula yields 

^o.{A;J,A,l)^^^co.iA,,A,,). (A-11) 

Hence the terms dA^J^/OAij are required. 
Consider the identity, 

A^'A,k = . (A-12) 
Taking the derivative with respect to Aa/3 yields. 

Making the substitutions. 



and 



yields 



dAjk 

dAaf3 



6jJkp (A-14) 



SiaSjb (A-15) 



Ar^'6,p + A,,-^ = 0. (A-16) 



Contracting with Aj^^ and using A^ii^Ai^^ — Si)^ gives 



which can be relabelled as 



dA 

~ ^ia A 3^ , (A- 17) 



-gf- = ■ (A-18) 
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Hence, equation [A-ll| becomes 

cov(A-^, A-J) = A-jA-^^A-lAi,'coY{A,„ A^i) • (A-19) 

This is the fundamental formula of this paper. 

An important case is when matrix elements Aij are uncorrelated. If each 
element of Aij has associated with it a variance [<J A^ij then 

cav {Aij, Aki) = [cTAl^jSikSji (no summation). (A-20) 

Note that a corrected version of equation ^ now can be calculated as the 
variance [cyi-i]^^ for each element of A'^j, 

Ml, ^ cov(A-^, A-J) = [A-']lAaA]UA^% • (A-21) 

Finally, the full covariance formula for uncorrelated Aij is given by 

coy{A-'^, a-,') = iA-lA-l)[aA]UA-,'^,') (A-22) 

where there is no sum inside the parentheses. For an n x n matrix there 
are covariances. From symmetry of the covariance operation one sees that 
there are ?t,^(?t,^ + 1)/2 independent covariance terms of which are diagonal, 
i.e. variances, and r?{v? — l)/2 are "off- diagonal" . 

B An Explicit Two Dimensional Example 

The two dimensional case is very instructive because it is simple enough to 
derive analytically. Consider 

where 

^ ^ ( - " ) "^-^^ 

and \ A \= ad — he Ys, the determinant of A. The matrix elements of A have 
uncertainties given by 
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The uncertainty on a = Aii, for example, can be calculated using the error 
propagation formula, 



The partial derivatives can be calculated from equation [B-2| as for example, 

^ - A ( ^ \ - — ( ^ 
(9a da \ad — be J da \ {ad — hcY ) 

For completeness the required derivatives are: 



a' . (B-5) 



da 2 n n /T-, \ 

and therefore 

2 4 2 I 2 2 2 I 2 <32 2 I <32 2 2 /t> ^\ 

This result can be compared directly with the results from equation A-21 



where a = P = 1, 

cov(Ar/, A^l) ^ = [A-%[crA]UA-% , (B-8) 

which on substitution for the elements of and cta yields equation P-7| . 

Similarly an explicit example of an "off-diagonal" covariance can be cal- 
culated, 

da dd r, dadd r, da dd ^ dadd ^ x 
where the partial derivatives with respect to beta are given by. 

Therefore the covariance of a and /5 is given by, 

cov(a, P) = a^pal + a^^5al + ap^al + P'^i5al . (B-11) 

This result is identical with equation |A-22| where a = /5 = l,a = l and 
6 = 2. Note, that in general, the covariance is not zero. 
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C Covariances of the Determinant 

Let A be an n X n matrix with elements Aij. The determinant is given by, 

det A = \A\ = €aia2...a„Alai^2a2 ■ ■ ■ ^na„ ; (^"1) 

where eo^Q,2,,,c„ = 1 for cychc permutations of ai, a factor of -1 is apphed each 
time any two indices are exchanged, hence eaia2...a„ is zero if any two indices 
are the same. This can be rewritten by defining new matrices by removing 
from A the elements of its i^^ row and j^^ column. The determinant of the 
remaining (n — 1) x (n — 1) matrix is called the minor of Aij. If the minor 
of Aij is denoted by M^- then the cofactor of Aij is given by, 

= i-iy^'M,,. (C-2) 

The matrix of all cof actors is C. The determinant of A can now be written 

as, 

det A = AijCij = AjiCji, (C-3) 

for i fixed. 

The error on det A is given by, 

2 9 det ^9 det A , ^ ^ , ,^ . 



From equation |C-3| one gets, 

a det A 



OA, 



Cab, (C-5) 



ab 



since z = a is fixed and Cij is independent of Aab for all j. Therefore, one 
obtains the very compact result, 

0-^41 = CabCaf3COv{Aab, A^p). (C-6) 
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